In this assignment, you will explore curvature on discrete surfaces. First you compute the Gaussian curvature and then the Mean curvature normal. Finally, you compute the shape operator for each vertex and use it to find both the mean curvature (again) and also the principal curvatures.
from pygel3d import hmesh, jupyter_display as jd
from numpy import zeros, dot
from numpy.linalg import norm
import numpy as np
from math import pi, acos, tan, acos, cos, sqrt
jd.set_export_mode(True)
# A couple of meshes are provided in the zip. Feel free to try more.
filename = "bunnygtest.obj"
# filename = "torus.obj"
is_bunny = ("bunny" in filename)
m = hmesh.obj_load(filename)
hmesh.triangulate(m)
jd.display(m)
Compute the Gaussian curvature at all vertices. Visualize the Gaussian curvature on the mesh and compute the integral gaussian curvature over the mesh. Print out the result of the integral Gaussian curvature. Try this both for the bunny and the torus. Note that when visualizing the curvature, some vertices are likely to have extreme values. You probably need to clamp the values to a range in order to be able to get a meaningful visualization.
def clip_worst(a,fraction=0.1):
a_sorted = np.sort(a)
n = len(a)-1
return np.clip(a, a_sorted[int(n*fraction)], a_sorted[int(n*(1.0-fraction))])
K = zeros((len(m.vertices()),))
pos = np.array(m.positions())
kint = 0
for v in m.vertices():
# Insert code here which computes the Gaussian curvature for each vertex,
# adds to the sum (kint), and assigns the value to K[v]
# Your code here
# ------------------
phi_sum = 0
area = 0
for fid in m.circulate_vertex(v, mode="f"):
# Sum all triangles areas
area += m.area(fid)
# Calculate phi with cosine rule
vid_con = np.array(m.circulate_face(fid, mode="v"))
v1,v2 = vid_con[vid_con != v]
a = np.linalg.norm(pos[v1] - pos[v2])
b = np.linalg.norm(pos[v] - pos[v1])
c = np.linalg.norm(pos[v] - pos[v2])
phi_sum += np.arccos((b*b + c*c - a*a)/(2*b*c))
K[v] = (2*np.pi - phi_sum) / area # Definition 8.3
kint += K[v]
print("Integral Gaussian curvature", kint)
if is_bunny:
K = clip_worst(K, fraction=0.1)
jd.display(m,data=K)
Integral Gaussian curvature 995864.8226074256
In this exericse, you should compute the mean curvature normal $\mathbf{H} = \text{H}\cdot \mathbf{n}$ for each vertex in the mesh. On the mesh, we plot the length of this vector which is the mean curvature.
def normalize(v):
return v / np.linalg.norm(v)
def vertex_areas(m):
pos = m.positions()
face_areas = [ m.area(f) for f in m.faces() ]
ora = zeros(len(pos))
for v in m.vertices():
edges = [normalize(pos[vn]-pos[v]) for vn in m.circulate_vertex(v, mode='v')]
areas = [face_areas[f] for f in m.circulate_vertex(v, mode='f')]
indices = range(len(edges))
angles = [acos(dot(edges[i],edges[(i+1)%len(edges)])) for i in indices]
ora[v] = sum([areas[i]*angles[i]/pi for i in indices])
return ora
def cot(theta):
return np.cos(theta)/np.sin(theta)
def spatial_elem(m, v, hid):
pos = m.positions()
# alpha is the angle opposite to hid
v1 = m.incident_vertex(hid)
v2 = m.incident_vertex(m.next_halfedge(hid))
a = np.linalg.norm(pos[v1] - pos[v2])
b = np.linalg.norm(pos[v] - pos[v1])
c = np.linalg.norm(pos[v] - pos[v2])
alpha = np.arccos((a*a + c*c - b*b)/(2*a*c))
# Beta is similar except for the opposite hid
hid_opp = m.opposite_halfedge(hid)
v2 = m.incident_vertex(m.next_halfedge(hid_opp))
a = np.linalg.norm(pos[v1] - pos[v2])
c = np.linalg.norm(pos[v] - pos[v2])
beta = np.arccos((a*a + c*c - b*b)/(2*a*c))
pi = pos[v]
pj = pos[v1]
return (cot(alpha) + cot(beta)) * (pi - pj)
def mean_curvature_normal(m):
N = m.no_allocated_vertices()
Hn = np.zeros((N,3))
# Insert code here which computes the Mean curvature normal H*n for each vertex
# using the contan formula.
ora = vertex_areas(m)
for v in m.vertices():
# Your code here
# ------------------
spatial_sum = np.zeros(3)
for hid in m.circulate_vertex(v, mode="h"):
spatial_sum += spatial_elem(m, v, hid)
Hn[v] = spatial_sum / (4.0 * ora[v])
return Hn
# Display Mean Curvature
Hn = mean_curvature_normal(m)
if is_bunny:
Hn = clip_worst(Hn, fraction = 0.2)
H = norm(Hn,axis=1)
print("Mean Curvature", H.shape)
jd.display(m, data=H)
Mean Curvature (3485,)
In this part, you compute the shape operator $\mathbf{S} = - \begin{bmatrix} a & b \\ b & c \end{bmatrix}$ and use it to find the principal curvatures $\text{kmin}$ and $\text{kmax}$ and the mean curvature $\text{H}$.
def shape_operator(m,v):
pos = np.array(m.positions())
# Insert code here which computes the shape operator matrix S that can be used to compute
# the principal curvatures kmin and kmax and the mean curvature H
# Your code here
# ------------------
# Create the a,b,n span
n = np.array(m.vertex_normal(v))
a_tilde_vec = np.random.rand(len(n)) # Create random array that is not parallel to the normal
while np.abs(n @ a_tilde_vec.T) == 1:
a_tilde_vec = np.random.rand(len(n))
a_vec = (a_tilde_vec - n * (a_tilde_vec @ n))
a_vec /= np.linalg.norm(a_vec)
b_vec = np.cross(n,a_vec)
# create U and F
v_neigh = list(m.circulate_vertex(v, mode="v"))
p_neigh_centered_at_v = pos[v_neigh] - pos[v]
T_T = np.vstack((a_vec, b_vec))
uv = T_T @ p_neigh_centered_at_v.T
U = np.vstack((uv[0]*uv[0]/2, uv[0]*uv[1], uv[1]*uv[1]/2)).T
F = p_neigh_centered_at_v @ n
a,b,c = (np.linalg.inv(U.T @ U) @ U.T) @ F
# Create shape operator
S = - np.array([[a,b],[b,c]])
S_3d = T_T.T @ S @ T_T
eig_vals = np.linalg.eigvals(S_3d)
index_nearest_zero = np.argmin(np.abs(eig_vals))
# Principal curvature
kmin, kmax = sorted(eig_vals[np.arange(3) != index_nearest_zero])
# Mean curvature
H = (kmin + kmax) / 2
return kmin, kmax, H
# Display Mean Curvature
H = zeros(len(m.vertices()))
for v in m.vertices():
kmin, kmax, mean_curvature = shape_operator(m,v)
H[v] = mean_curvature
H = clip_worst(H, fraction=0.1)
print("Mean Curvature")
jd.display(m,data=H)
Mean Curvature
# Display Principal Curvatures
Hn = mean_curvature_normal(m)
min_curvature = zeros(len(m.vertices()))
max_curvature = zeros(len(m.vertices()))
for v in m.vertices():
kmin, kmax, mean_curvature = shape_operator(m,v)
min_curvature[v] = kmin
max_curvature[v] = kmax
min_curvature = clip_worst(min_curvature, fraction=0.1)
max_curvature = clip_worst(max_curvature, fraction=0.1)
print("Principal curvatures")
jd.display(m,data=min_curvature)
jd.display(m,data=max_curvature)
Principal curvatures
Answer: The positive curvatures are colored in green and the negative curvatures are colored in pink.
Answer: Yes, mean curvature would take the mean of the curvatre to all directions from a point. We would probably not see so obviously the sharp changes between the inner and the outer parts.
Answer: The value gives us information about the shape of the object. This is basically the sum of the total curvature of the object. The bunny is all just curving "out" so the curvatures are positive in most points so the value is big and positive. The Torus has hole and many points have negative curvatures (also many points with positive). So the value for the torus is much closer to zero.
Answer: The outside of the Torus is a positive curvature and the inside of the torus is a negative curvature. The mean curvature is always positive is the average of both the positive and the negative, but the outside part influences the mean more and therefore the mean is always positive for the torus.
Answer: On a torus, the maximum curvature stays the same because it's like the thickness of the torus, that doesn't change. The minimum curvature changes a lot because it is like going from the inside of the hole to the outside, and that shape changes a lot.